DATA_SECTION
  init_int Nobs_tot                    //Total number of observations
  init_vector Kw(1,8)                  //Value of Water absorption for each wavelength
  init_matrix Data(1,Nobs_tot,1,6)     // Data arranged in column such as :Kd, Chl, Acdm, Bbp, lambda, convert_lambda
  vector Kd(1,Nobs_tot)                // Kd data as a vector
  vector Chl(1,Nobs_tot)               // Chlorophyll data as a vector
  vector Acdm(1,Nobs_tot)              // Acdm data as a vector
  vector Bbp(1,Nobs_tot)               // Bbp data as a vector
  vector lambda(1,Nobs_tot)            // Lambda as vector
  vector convert_lambda(1,Nobs_tot)    // Lambda as an indice 1=320, 2=340, 3=380, 4=412, 5=443, 6=490, 7=510, 8=555

PARAMETER_SECTION
  init_vector logSigma(1,8)
  vector sigma(1,8)
  init_number s0
  init_number s1(-1)
  init_number s2(-1)
  init_vector loga1(1,8)
  init_bounded_vector a2(1,8,0,1.01)
  init_number eta
  vector pred_Kd(1,Nobs_tot)	       		//one vector of the predicted value of Kd
  sdreport_vector a11(1,8)
  objective_function_value f

PRELIMINARY_CALCS_SECTION
  Kd=column(Data,1);
  Chl=column(Data,2);
  Acdm=column(Data,3);
  Bbp=column(Data,4);
  lambda=column(Data,5);  
  convert_lambda=column(Data,6);

PROCEDURE_SECTION
  f=0.0;
  sigma=exp(logSigma);
  a11=exp(loga1);
  for (int i=1;i<=Nobs_tot;i++)
  {
  pred_Kd(i)=Kw(convert_lambda(i))+Acdm(i)*mfexp(-(s0+s1*Acdm(i)+s2*pow(Acdm(i),2))*(lambda(i)-443.0))+a11(convert_lambda(i))*pow(Chl(i),a2(convert_lambda(i)))+Bbp(i)*pow((lambda(i)/443.0),eta);
  f+=0.5*(Nobs_tot*log(2*M_PI*square(sigma(convert_lambda(i))))+square((log(pred_Kd(i))-log(Kd(i)))/sigma(convert_lambda(i))));
  }
//  f=0.5*(Nobs_tot*log(2*M_PI*square(sigma))+norm2((log(pred_Kd)-log(Kd))/sigma));

RUNTIME_SECTION
  maximum_function_evaluations 200000
  convergence_criteria 1.e-12

TOP_OF_MAIN_SECTION
  arrmblsize = 100000000;
  gradient_structure::set_MAX_NVAR_OFFSET(30000);

